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We study Mode I fracture in a viscoelastic lattice model with a nonlinear force law, with a 
focus on the velocity and linear stability of the steady-state propagating solution. This study is 
' a continuation both of the study of the piece-wise linear model in Mode I, and the study of more 

general nonlinear force laws in Mode III fracture. At small driving, there is a strong dependency of 
■ the velocity curve on the dissipation and a strong sensitivity to the smoothness of the force law at 

large dissipation. At large driving we calculate, via a linear stability analysis, the critical velocity 
for the onset of instability as a function of the smoothness, the dissipation and the ratio of lattice 
spacing to critical extension. This critical velocity is seen to be very sensitive to these parameters. 
We confirm our calculations via direct numerical simulations of the initial value problem. 



(N 



PACS numbers: 62.20.Mk,46.50.+a 

CZ2 



O 



I. INTRODUCTION 



The problem of extensional (Mode I) fracture has received increasing attention in the physics community in the last 
decade J lj. The experiments in amorphous materials done by Fineberg, et al. Q and in ordered materials by Hauch, 
ct al. |3|, showing interesting dynamical behavior for large velocity cracks have been at the center of this growing 
interest. From a theoretical point of view, the singularities at the crack tip associated with continuum treatments of 
the crack make the problem challenging. The presence of these singularities necessitates a treatment of the crack tip 
at the microscopic level, the singularities being regularized by the small scale dynamics. One line of attack on this 
problem, initiated by Slepyan [Q, has been through the study of lattice models of cracks. These lattice models are 
simpler than full atomistic simulations Q in that the connectivity of the atoms is specified from the beginning, and 
fSJ \ so dislocations are excluded. 

■ Much progress has been made in understanding steady-state propagation of cracks in lattice systems using an 
ideally-brittle piece- wise linear force law [Q, g, 0, [| 11, Q. Here the particles interact with Hookean springs, 



which break at some critical extension, after which they exert no force. With these piece-wise-linear interactions, the 
model admits an analytic solution via the Wiener-Hopf technique. This solution has been carried out both for Mode 
III and Mode I cracks, for both finite width and infinitely wide systems, with and without dissipation (Stokes or Kelvin 
type). There have also been some recent results on mode II cracks, with possible relevance for an understanding of 
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The piece- wise linear model, despite its analytic simplicity, exhibits some undesirable features For small 

dissipation, the solutions are inconsistent at small velocities, in that bonds on the crack surface which are assumed 
to crack at time some specific time in fact are seen to reach the critical displacement and hence to crack earlier. At 
large velocities, the analytic solutions are again inconsistent, this time due to the breaking of additional bonds off 
the nominal crack surface. Beyond this point, the analytic methods are unable to tell us anything about the true 
dynamics of the system. An additional limitation of the piece- wise-linear force law is that it greatly complicates the 
£j \ task of constructing a linear stability treatment of the steady-state crack. This is due of course to the discontinuous 
r* ■ nature of the force law. 

For the case of Mode III fracture, Kessler and Levine jl3| studied a lattice model with a family of continuous 
non-linear force laws labeled by a tunable smoothness parameter a. Such continuous forces should presumably more 
closely model the actual physical situation. With this lattice model, they investigated both arrested cracks |l4| ] 
and propagating cracks 13 . Consistent steady-state solutions were found at all driving; the previously discovered 
inconsistent solutions at large driving of the piece- wise linear limit were converted to consistent but linearly unstable 
solutions for the continuous force law. 

In this paper we will study Mode I fracture with a suitable adaptation to vector forces of this class of continuous force 
law. We shall study steady-state solutions of the model both for small and large driving, and their dependence on the 
dissipation, the critical bond extension (relative to the lattice constant) and smoothness. Our results will be compared 
those of both the piece- wise linear limit Mode I calculations jl(J] and the continuous force Mode III analysis jl^] . We 
shall also study the instabilities of the theory at large driving, finding a number of transitions between dominant 
modes as a function of dissipation; these findings will be confirmed by direct numerical simulations. We will also 
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briefly consider the nature of the dynamics, post-instability. Finally, we shall draw some conclusions and discuss 
directions for further research. 



II. MODEL AND GENERAL METHODOLOGY 



In this paper we study extensional (Mode I) cracks in a nonlinear, two dimension hexagonal lattice model. Mass 
points located at the lattice sites are coupled to their six nearest-neighbors by nonlinear, viscoelastic central forces. 
The force that particle 2 exerts on particle 1 is taken to be: 

r , s 1 + tanh[a(e - ri, 2 )] , . 

fl ' 2 = (ri ' 2 a) l + tanh(a) ^ (1) 

where a is the lattice scale, r 12 = X\ — X2, ri.2 — is the distance between the atoms and f 12 = The 

threshold e is the distance between the atoms at which the spring can be considered to break, a is a parameter 
that determines the smoothness of the potential. In the limit a — > 00 the force becomes perfectly brittle, dropping 
immediately to as ri.2 exceeds e. Decreasing a smoothes the force, so that it decays away over a distance 1/a. 
Now the bonds never crack totally, there is always some (exponentially small) force, even at large ri,2- Note that, 
compared with the similar force law studied at Mode III [13], this model has an additional parameter, namely a/e the 
ratio of the lattice constant to the breaking extension; this parameter governs the extent to which the force direction 
changes as the displacement grows. Only in the joint limit a — > 00, a/s — > 00 do we recover the piece- wise linear force 
law studied in Mode I by Kulamekhtova, Saraikin and Slepyan Marder and Liu (l6) and Pechenik, Kessler and 
Levine |K^. In the following we choose our length scale such that s = 1. 

In addition to the purely conservative force, we introduce a Kelvin-type viscosity with a viscosity parameter r\ with 
the force law f§ |, [h], |l|: 

§1,1 = #1,2(^1,2 • n,2)n,2 (2) 

where «i 2 = V\ — i/ 2 and fci j2 is an effective spring constant: 

ki,2 = /i,2/(ri,2 - a). (3) 

These forces define our model. We can now write the equation of motion for the displacement u of a mass point away 
from its lattice position: 

d 2 u{x) ^-^ , -» _ 

dt 2 = U 3,2' + 93,3') ■ (4) 

x' Gnn 

We work on a hexagonal lattice with 2N + 2 rows in the y direction, indexed by j = —N . . . N + 1. The rows are 
separated by a distance \/3/2 • a, so that yj = (2j — l)a\/3/4. We apply a constant displacement u = ±Ay to the 
top- and bottom-most rows. The (metastable) uncracked state is that of uniform strain 

^(x, yj ) = ^±Ay. (5) 

In the equilibrium cracked state, the upper half of the rows have u — +Ay, and the lower half u = —Ay. 

Initially, we will be interested in the case of steady state cracks where the displacement has the Slepyan form: 
u(t,x,y) — u(t — x/v,y) where x,y label the position of a particle on the lattice, and v is the crack propagation 
speed ||, This ansatz reduces the problem to one of solving for the time development of the displacement for 
just 2N particles, i.e. one particle on each of the unconstrained rows. In addition, we focus on symmetric cracks so 
we need consider just one side of the lattice, imposing the symmetry condition: u v (t,y) = —u Vt {t — a/(2v),—y) and 
Ux(t, y) =u x (t- a/(2v),-y). 

As opposed to the previous analyses done in the piece- wise linear limit (a, a — * 00) fl(i|| , here with finite a, a we 
must resort to a numerical procedure. Specifically, we discretize time using a small but finite dt, and further limit 
— T ^ t ^ T. This reduces the problem to a set of 2N(2T + 1) coupled equations. Since the equations of motion are 
dependent only on nearest neighbors, the system has a banded structure. The equations of motions are nonlinear so we 
use Newton's method to solve the problem, using standard subroutines for solving banded matrices to find the Newton 
update. Now, for each given v of the crack, we want to find the driving displacement A. We have 2N(2T + 1) + 1 
unknowns (the N vector functions for the entire time range plus A); this is one more than the number of equations, 
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FIG. 1: (a-c) Dependence of v/cr on A/A G for a = 5, 15 and 100, for r] = 0.25, 1 and 3. N = 3, dt = 0.1 and a = 4. In (c) 
we present data also for a — 8. 



reflecting the time-translation invariance of the sought-after solution. Thus we add an additional equation fixing 
Uy_/y^(0) — const to lift the degeneracy. Let us note that it can be shown that the discretized problem has one 

more linear mode that diverges as x — > — oo than diverges at x — > oo. To treat this, we enforce the equation of motion 
for — T — dt ^ t ^ T — dt whereas the variable run from — T ^ t ^ T ||. The system then has 2N(nt + 3) — 2 bands 
below the diagonal and 2N(nt, — 1) + 3 upper bands, where rif, is a even number with v = a/(ribdt). 

The only remaining technical issue is that the system as presented does not have a completely banded structure 
due to the fact that many equations depend on A. This problem can be circumvented using the algorithm of Kessler 
and Levine (see the appendix of Ref. p3[). 



III. THE SMALL VELOCITY REGIME 



It is reasonable to expect that smoothing out the potential can have a large effect in the small velocity regime. In 
particular, the fact that solutions in the piece-wise linear limit become inconsistent at small enough velocity means 
that making a finite must make a significant difference. For mode III ([j~3|]), decreasing a eliminated the oscillations 
in the v(A) curve, in accord with this reasoning; but, the oscillations were already not present at large enough r\ and 
thus changing a has almost no effect on those steady-state curves. 

Here for Mode I we find a somewhat different picture. First we present in Fig. |l|(a-c) a graph of v versus A for 
the steady state solution, v is normalized with respect to the Rayleigh surface wave speed, the limiting crack speed 

(in the N — > oo limit) as the elastic field has to propagate along the crack surface: Cr = ^ 3 ~^ a ~ 0.563 a. A is 
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normalized by the Griffith value, Aq, the value for which fracture is energetically allowed. This value is a function 
of a, a and of course N. Data is presented for N = 3 because computing the very low velocity solutions is very hard 
at large N, due to the bandwidth which scale as N/v. We do not expect the picture at larger N to be qualitatively 
different @. 

At small 77 the curves have oscillations for all a, unlike the Mode III case, where oscillations were present only for 
large a. At large a, in both cases the graphs are qualitatively similar to those obtained directly in the piece- wise linear 
limit (PLL), though the large a curves do not actually converge to those curves. This is due to the aforementioned 
fact that the solutions of the PLL are actually inconsistent for small velocity whereas the solutions of the finite a 
model are always well-defined @, 1|. As we increase rj the oscillations disappear and the curves become smoother. 
This behavior, as can be seen in the 77 = 1 data (Fig. 0(b)), happens for all a, but is especially pronounced for the 
smaller a's as expected. These trends are also seen in the Mode III case |l|. As we increase r\ to still higher values, 
for example r\ — 3, (Fig. [j](c)), a maximal velocity can be seen at small a. The existence of a maximal velocity 
in steady-state cracks is a general feature, and will be discussed at length in the next section on the large velocity 
regime. Here we just note that for large rj and small a, the maximal velocity can be quite low, e.g. v max ~ 0.25C# 
at a = 5, a = 4, rj = 3, and decreases for decreasing a. This behavior has no parallel in the Mode III case. At small 
A, decreasing a increases the velocity, but for large A the situation is reversed, with the curves with the smaller a 
having lower velocity. This is a result of the fact that in our nonlinear potential the bonds never really break and 
continue to impede the crack's progress. On the other hand, at small A, small a weakens the lattice trapping effect, 
resulting in larger velocity than for larger a. 

We still have to address the question of why the curves oscillates at low a for small rj, in contradistinction to 
the Mode III case It is interesting to consider (see Fig. |J) the bond extension and the bond force for a 

representative small velocity case, v — 0.148Cr. In the piece-wise linear case at small rj for both Mode I and Mode 
III the underdamping of the backward running waves leads to pre-cracking. This means that the bond extension rises 
above unity (e = 1) before t = and then comes back to unity before cracking once and for all j?], ||. For finite a, 
while there cannot be any pre-cracking (we don't have to postulate in advance when a bond will break), there is a 
remnant of this tendency. After cracking, the bond extension first rises and then falls to a value close to the critical 
extension before rising again and then permanently leaving the region close to critical extension. This force oscillation 
is ultimately what is responsible for the oscillations in the v (A) plot. This situation is very different from that in 
Mode III, where for finite a, the bond extension increases monotonically through the critical region. 

We checked also the effect of the lattice scale constant a on the curves for the small velocity regime at 77 = 1, 
a = 15. We present the data on Fig. |. We see that there is no dramatic effect of the lattice scale constant in the 
small velocity regime. There is a backward branch for all a at about the same velocity, with the effect accentuated at 
large a. 



IV. LARGE VELOCITY REGIME 



In the large velocity regime, it is again clear that smoothing the potential must make significant modifications to 
the steady state solutions, as it is known that in the piece- wise linear model the solution is inconsistent due to the 
breaking of bonds off the assumed crack surface. This inconsistency was demonstrated both for Mode III in a square 
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FIG. 4: (a) v/cr versus A/ Ac for N — 10, r\ — 0.25 and a — 4. Data is presented for the cases a = 5, 7, 15, 30 and 100. (b) 
f/CR versus A/Ag for AT = 10, rj = 0.25 and a — 15. Data is presented for the cases a = 4, 10 and 75. 



lattice and Mode I in a hexagonal lattice, the case treated herein ffl, [H], 16 1. In Mode III, the inconsistency in the 
PLL is converted to a linear instability for the finite a model. 

We present in Fig. [|(a) large velocity steady-state curves for the case N = 10 for several values of a. The effect of 
changing a is rather dramatic. For a — 30, the velocity grows smoothly with A. For a = 100 the curve still grows 
but a little bit less smoothly with some small kinks appearing. There is no sign of a maximal velocity, at least for the 
range of A investigated; the maximal velocity, if it exists, is close to the Rayleigh wave speed. As we go lower in a, 
at a = 15, and more significantly at a = 7 and a — 5, the curves lie significantly below those for the larger a. We see 
here a maximal velocity far below the Rayleigh wave speed. These changes are a direct sign of the involvement of the 
additional bonds in the crack progression; in particular, the creation of a local maximum of the velocity as a function 
of A is due to the fact that this involvement grows with A and can overcompensate for the natural tendency of the 
crack to speed up as the applied stress increases. 

In the previous figure, the lattice scale a was kept fixed. If we increase a, the maximal velocity decreases and is no 
longer near the Rayleigh wave speed even at larger a. This is presented in Fig. |^(b), where we show the curves for 
a = 15 for different a. Thus the maximum velocity is a strongly varying function of both a and a. It will turn out 
from the linear stability analysis to be presented in the next section that there exists an instability which typically 
sets in slightly after the maximal velocity point; this means that the instability point will also be a strongly varying 
function of the parameters. In the experiments, it should be noted, the critical velocity for the onset of instabilities 
appears to be material dependent [Q . Thus, this strong dependence is an encouraging sign. Also, for large a and a, 
the maximal velocity is very close to the velocity at which the PLL becomes inconsistent, as expected. 
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V. STABILITY ANALYSIS 

Given the steady-state solution, we can proceed to investigate its linear stability. We will focus here on the instability 
in the high velocity regime; it is clear that in the low velocity regime where the v(A) curve is backward, the solutions 
are unstable. We take the form: 

Ui{t) =uf\T) + e^5ui{r) (6) 

where r is the traveling wave coordinate r = t — x/v, and u[°\t) is the steady state displacements for row i that 
we have found in the last chapter. We assume < 1 and expand the equations of motion to linear order. Note 
that e UJX / a — e UJV ( t ~ T )/ a - thus, stability requires Re(w) < 0, and stable perturbations decay ahead of the crack in the 
moving frame of reference. 

Substituting Eq. (jfj) into the equations of motion, Eq. (^) we get the linear equations: 

^Suiin) = £ e-*^ I fui ■ — (|— (./>.,• + //V .,• ) (7) 

1 jtnn \ uu j \>3) / 

where / and g are the forces defined in chapter II and Tj = t — Xi/v. 

To find the possible values of u> we proceed as follows: First we impose the boundary conditions 5u x ^ y = 0. for the 
constrained rows i — —N, N + 1. Also 5u x _ y = for any x's that correspond to times before or after the limiting 
values of r that were selected to build the steady-state solution. Second, we add a normalization condition by setting 
y ■ 5u\(0) = 1.; there are now as many equations as variables if one includes ui among the latter. Then, as in the 
steady state solution, we need to solve a set of (now complex) linear equations with an almost banded structure, where 
now uj is the variable that destroys the banded structure. To proceed, we choose some (complex) value of u), and 
temporarily ignore the actual equation of motion for the y component of Sui(0), replacing it by the aforementioned 
normalization condition. We then solve the banded system of linear equations by standard techniques. We still have 
left the equation of motion we have dropped, which serves as a complex mismatch function; the zeroes of this mismatch 
function determine the eigenvalues of u>. We are of course interested in the eigenvalue with largest real part, which is 
the dominant mode. 

When we solve the system of equations for the linear stability analysis, we do not impose any symmetry on the 
problem; i.e. we use the complete lattice encompassing both sides of the crack. Let us note that any time we find a 
complex root of to, with eigenfunctions Su, there is always an equivalent solution generated by the transformation: 

Imw — > 2n — Imw 

Su x ^y — > 5u X:V (even rows) (8) 
Su x>y — > — Su x _ y (odd rows) 

By even and odd we mean the serial number of the row, for example row is the first even row. This alternation 
is due to the shift of a/2 for the x- values in neighboring rows. 

To test our stability analysis we check whether our linear equations reproduce the time translation mode that should 
be always be present with eigenvalue lo = 0, and eigenfunction 

^(r) = ^uf ( T )/j-(y.^\0)). (9) 

Given our discrete St, we found the eigenvalue of the translation mode to lie not exactly at zero, but still quite small; 
e.g., typically uj ~ 10~ 4 for St ~ 0.1. Note that compared to the case of Mode III, the matrix here is larger as we have 
two components of the displacement field Furthermore, the critical velocity here is much smaller than the Mode III 
case, so n& is larger. It was therefore difficult to work at N = 10 as we did for the steady state solution. The stability 
matrix is complex as opposed to real, and we do not impose any symmetry on the problem; hence the matrix is now 
four times bigger than in the steady state solution case. We chose to work at N = 4 to reduce memory consumption. 
N = 4 is of course not good enough to give us quantitative results for a macroscopic system, but we expect the 
qualitative picture to remain the same. 

We focus first on the case of a = 15, a = 75, and 77 = 0.25. We chose a large value of a because the curves are much 
smoother than for small a. In addition, this choice allows for a more meaningful comparison to the results for the 
critical velocity for the PLL ETJ In Fig ||(a) we present the steady-state v(A) curve, along with the corresponding 
curve for N — 10. The overall similarity of the two curves, with an overall shift in velocity of about 5%, lends support 
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FIG. 5: (a) The v(A) curve for the case of TV = 4, <r) = 0.25, a = 15 and a = 75. We also show the N = 10 curve with the same 
parameters to demonstrate the effect of changing TV. (b-c) Eigenvalues for the stability problem, with TV = 4, r\ — 0.25, a — 15 
and a = 75. The calculation was done with 2T + 1 = 1000 and n;, = 20. Re(oj) is shown in (b), Im(u) in (c). 



to our expectation that the TV = 4 system is qualitatively similar to those with larger TV. In Fig. ||(b) we exhibit 
the real parts of the two most important eigenvalues as a function of A, and in Fig. ||(c) the imaginary parts of the 
eigenvalues. We can see by looking at the steady state solution curve and the linear stability curve that the point 
of instability lies slightly after the maximum in the v(A) curve. This is reasonable because both the instability and 
the maximum are signs of additional bond breaking, and so should occur close together. In the next chapter we will 
verify this behavior via direct simulation. The second interesting point is that in all our solutions there is a stable 
mode that almost does not change with A, and there is a quickly- varying mode that takes its place as the dominant 
mode and crosses over to positive Re(cj). 

In Fig. ^ we show a typical eigenfunction (for the quickly-varying mode) from the linear stability analysis. Fig. 
||(a) presents the real part of the y component of the eigenfunction for the two rows on either side of the crack surface 
(y = ±\/3/4) as a function of the traveling wave coordinate r. We see that in this case the most unstable eigenfunction 
is antisymmetric about the crack surface. The eigenfunction decays slowly downstream and very rapidly upstream of 
the crack tip. In Fig. ^(b) we show the real part of the x component for the two rows about the midline. Here, both 
upstream and downstream of the crack tip, this component of the eigenfunction decays very slowly. Moreover, it is 
here symmetric about the midline; the two curves overlap. We note in passing that there is a large class of possible 
symmetries for the eigenf unctions. In addition to the eigenfunction with antisymmetric y component and symmetric 
x component, there also exists eigenf unctions with symmetric y component and antisymmetric x component. The 
symmetry of the imaginary part of the eigenfunctions is also variable, and not connected in any special way to the 
symmetry of the real part. Thus, while imposing a specific symmetry reduces the size of the matrix to be solved, it 
means that the calculation has to be repeated for each type of possible symmetry. 

We next investigated a higher value of T), T] = 0.75. The steady-state curve, as we see in Fig. 0(a), rises to a 
maximum and then it comes back around in a lower branch. In the linear stability analysis for this 77, Fig. |?](b), we 
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FIG. 6: Eigenfunction for stability problem, with N = 4, r\ — 0.25, a — 15 and a = 75: (a) y component, Re y ■ 5uo,ij (b) 1 
component, Rey • 5uo,i. Both figures were done with 2T + 1 = 1000, nt, = 20 and dt = 0.13. 

see that the point of instability is exactly at the A-maximum; all the points on the backward branch are unstable. This 
is of course a generic feature of propagating systems. At the A maximum, the velocity of the solution is indeterminate 
to linear order, implying a zero mode at this point. This zero mode represents the crossing over into instability of 
some mode of the system. When we decrease 77 the backward branch comes back less rapidly, as shown in Fig. @ for 
•q = 0.5. Clearly at some smaller 77 this turnaround vanishes altogether, and we recover the situation at 77 = 0.25, 
where no maximum A is seen. The linear stability analysis for the r\ — 0.5 reveals that the maximal A is no longer 
the location of the first instability; there is another mode which goes unstable before this point. 

We present in Fig. 51(a) a curve of the critical velocity for the onset of the instability as a function of 77. Given the 
correlation we found between the critical velocity for instability and the point of maximal A for sufficiently large 77, 
we extended the graph to higher 77 by calculating the velocity at the A maximum. The critical velocity for small 77 
increases with 77, as in the piece-wise linear model (where here the critical velocity denotes the onset of inconsistency) ; 
even the values are similar. As we increase 77 further, the critical velocity decreases, again just like the piece- wise linear 
results, due there to the breaking of a horizontal bond fIo|| . We further see that the critical velocity extrapolates to a 
finite value at 77 = 0. Thus, there is stable propagation below this critical velocity at zero dissipation in contradiction 
to the claim of Pla, et al. |17| ■ We believe that the instabilities they saw in their simulations were finite amplitude 
effects engendered by the initial conditions. Indeed, we have generated stable propagating solutions at 77 = via 
direct simulation (see next section). 

Note that in the Mode I case as opposed to the Mode III case Im(tj) at the point of instability is not close 
to 7r. Therefore, there is no parallel to the period doubled nonlinear state found in that system. The critical Im(cj) 
here varies strongly with 77. We present in Fig. ||(b) the hxi(u>) at the point of instability as a function of 77. We see 
two transitions between dominant modes; somewhere in the range 0; .78 < 77 < 0.83, the critical Im(aj) becomes zero, 
as it must if the transition is directly connected to the saddle-node behavior of the steady-state curve. At lower 77, 
e.g., 77 = 0.5, it was about 2.0, as mentioned above. Note that this transition is sharp; for example at 77 = 0.83 we 
still have the mode with lm(u) ~ 2 but the Im(w) = mode goes unstable first. We also have a transition between 
dominant modes at 77 ~ 0.4, the modes having very different Im(w). At higher 77 we checked the critical Im(w) by the 
direct simulation (again, see the next section) and found that it remains zero. 

VI. COMPARISON TO SIMULATION 

To check our results, both for the steady state solution and for the linear stability analysis we have implemented 
a direct numerical simulation of the initial value problem. These simulations are also useful for investigating the 
dynamics after the point of instability. Our system is a hexagonal lattice with 2N + 1 rows and with a finite number, 
L, of mass points in each row. The top and bottom rows are constrained to have a displacement ±Ay. The initial 
conditions for the displacement and velocity for the lattice were constructed from the steady state solution that 
we already have found, with the initial crack tip placed L/3 from the left edge. We solve the equations of motion 
by an Euler scheme, taking Vt+dt = Vt + ^""^(/ + (f)dt and x t +dt = s?t + vt+dtdt. We measured the velocity of the 

nn 

crack by monitoring the bond extension. The crack was considered to have translated forward when both bonds 
connecting a given mass point to its neighbors across the midline exceeded the critical extension, e = 1, and the 
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FIG. 7: (a) The v(A) curve for the case of N — 4, rj = 0.75, a = 15 and a = 75. (b) Re(u) of eigenvalue for stability problem, 
with N — 4, rj = 0.75, a = 15 and a = 75. The calculation was done with 2T + 1 = 1000 and m = 20. 
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FIG. 8: The v(A) curve for the case of N = 4, rj = 0.5, a = 15 and a = 75. 

velocity determined from the time At between such translation events, v = a/ At. This criterion of velocity is of 
course somewhat arbitrary because cracking in our smooth potential is a reversible process. Note that to produce 
a true steady state crack we need to fix a/v so that it is an integer multiple of dt. Otherwise we would introduce 
oscillations in At due to the incommensurability. Via this procedure, we obtained excellent agreement (with 4 digit 
accuracy) between the simulated and calculated results for the v(A) relationship. As we increase A, we can follow 
the v(A) curve and we obtain from the simulation all the points of the curve including those around the ^-maximum. 
When we get exactly to the point of instability as predicted from the linear stability analysis, additional bond breaking 
occurs, and the velocity in the simulation drops down below the steady-state curve. 

The numerical simulation can provide a strong check on the details of the stability analysis. We can measure 
the eigenfrequency by perturbing one mass point on the lattice in the vicinity of the crack tip, and measuring the 
subsequent time development of the perturbation. We track the bond extension of the breaking bond at the moment 
of cracking. Because of the finiteness of dt this is always somewhat larger than e. This value oscillates, decaying in 
the stable regime, and growing in the unstable regime. The data can then be fit to a form Ae^ and ui compared 
to that produced by the linear stability analysis. In Fig. |l(](a) we can see the results for the crack tip and if we 
perform the fitting, the eigenvalue found here lies quite close to the eigenvalue found in the linear stability analysis, 
uj = —0.0243 + 0.5813z. If we instead look at a point in the unstable regime, this oscillations look like those presented 
in Fig. |l^(b). The oscillating part in the simulation fits well to the imaginary part of the linear stability eigenvalue. 
We note that throughout the unstable regime we can detect the presence of the square of the unstable mode, giving 
rise to an additional non-linear mode whose real part is two times the value of the real part of the eigenvalue. 

When we increase A past the point of instability, the actual crack dynamics of the problem become increasingly 
complex. Direct numerical simulation an still be used to study crack development. We performed such a study using 
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FIG. 9: (a) The critical velocity from the linear stability analysis, normalized by the Rayleigh wave speed, as a function of r\. 
Data is presented for a — 15, a = 75 and TV = 4. (b) Im(o;) of the eigenvalue at the critical velocity as a function of r\. 
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FIG. 10: (a) Bond extension at breaking as a function of position x along the crack surface. Here N = 4, r) = 0.25, a — 15 
and a = 75. A/A G = 1.6305, dt = 0.125. The fitted curve is 1.021 + COle -0 025 * sin(0.565a; + 3.5). (b) Bond extension at 
breaking as a function of position x after the point of instability for N = 4, r) = 0.25, a = 15 and a — 75. 



N = 20. We see in Fig. [Ll] three examples of what can happen in Mode I cracking after the point of instability, for 
increasing values of A, all for the case of small dissipation. For the smallest A, the crack tip meanders up and down 
off the midline. For larger A, the crack starts to bifurcate (Fig. |ll|(b)), and at the highest A shown bifurcates many- 
times and creates a very complex cracking picture (Fig. pd|( c)) . In all three of these small 77 cases, the crack arrests 
eventually. At large 77, such as rj = 1.2 and higher (see Fig. |lj) the picture is very different. Just above the point of 
instability the crack bifurcates, with the two daughter tips propagating close to the edges of the system. These kinds 
of behavior appear to be very different from that seen in experiment, where there is a main crack propagating more 
or less straight in addition to daughter cracks on either side. 

In the Mode III case, a crucial point that emerged from the detailed linear stability study was that the additional 
cracking bonds were always created behind the crack tip. There, the instability did not divert the main crack, but 
merely slowed it down temporarily as energy was diverted to create the sidebranches |l3|, [l8| . Here in Mode I cracking, 
however, after the point of instability the crack tip itself changes direction and the crack either bifurcates or deviates 
from the center. 



VII. SUMMARY AND DISCUSSION 



We have herein studied the steady-state crack, calculating the v(A) curves both for small and large driving, and 
for small and large dissipation. For small driving with small dissipation, we obtained oscillating curves for all finite 
a, indicating a substantial velocity gap at small velocity. At large dissipation, the curves are smooth, and vary 
significantly with a. At large driving, the maximal velocity exhibits a strong dependence on the smoothness and the 
value of a/e. For large 77, there is also a maximal driving as the curve bends back forming a second branch, and 




FIG. 11: (a) The lattice showing the broken bonds for N = 20, 77 = 0.1, a = 15, a = 75. A/A G = 1.875. (b) A/A G = 2.1875. 
(c) A/A G = 3.125. 
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FIG. 12: The lattice showing the broken bonds for N = 20, Tj = 1.2, a = 15, a = 75. A/A G = 2.578. 

beyond which no steady-state solution exists. For all 77, the steady-state solution undergoes a linear instability. For 
large 77, this occurs at the maximal driving, and is a real instability. For small 77, the instability set in slightly after 
the maximal velocity, and is a Hopf bifurcation. The velocity at the onset of instability initially increases with 77, 
reaches a maximum at 77 ~ 2 and thereafter decreases for larger 77. We note that the strong variation of the velocity 
at the onset of instability with the microscopic parameters does not accord with the naive expectations based on the 
Yoffe continuum calculation [[l9[ of a change in the direction of maximal stress at some universal velocity, independent 
of the microscopic dynamics. It also is not accord with the ansatz of Eshelby based on energy considerations [p0| . 
We also compared our steady state solution and linear stability analysis to direct numerical simulations, obtaining 
excellent agreement. 

We plan to extend our work to the case of biaxial loading, where the material in strained in both the x and y 
direction. This phenomena describes for example a popping of a balloon as described in Ref. pi] ]. We also plan 
to investigate further the post-instability dynamics. Finally, we hope to use the recently introduced p2| continuum 
rcgularization of tip-dynamics (based on the phase-field method) to unravel the role that the lattice structure has in 
determining the form of the allowed instabilities. 
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